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Abstract 



. Using a hydrodyamic model of granular flows, we present very 

p^ . long time simulations of a granular fluid in two dimensions without 

C^^ i gravity and with periodic boundary conditions in a square domain. 

/Cj I Depending upon the values of the viscosity, thermal conductivity and 

^^ ■ dissipation, we find for intermediate times a metastable clustering 

^J , state. For longer times the system is attracted to either a shear band 

or a vortex state. Our results are in general agreement with molecular 



C^ ■ dynamics simulations. 

"§ ■ 1 Introduction 

o 

^ ' Granular materials are an economically important and physically interest- 

ing systemPP. These materials are dissipative in nature and have complex 
^ ■ behaviour and exhibit very interesting phenomena such as clustering P| and 

Maxwell Demon effects^]. The fact that the behaviour of these materials 
is in many ways similar to conventional fluids was noted very early and a 
hydrodynamic approach to them was developed heuristically by Haff |3] and 
Jenkins and Savage [Sj. More recently a similar hydrodynamic approach has 
been developed by considering a Chapman-Enskog expansion of the Boltz- 
mann equation for a granular system by Brey, Dufty, Kim and Santos jHl El- 
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In this paper, we study a modification of tlie Haff equations by Hill and 
MazenkofH] which enabled the use of the hydrodynamic equations to be ex- 
tended to systems with clustering instabilities. This model has the advantage 
of avoiding inelastic collapse and of being able to be simulated efficiently on a 
computer. Further, it was also shown by Hill and MazenkojH] that these hy- 
drodynamic equations form a well-defined system of equations under a range 
of circumstances. It therefore also suggests the possibility that the nonlinear 
hydrodynamical approach is applicable outside the conventional regime of low 
densities. We investigate this model in detail and show that it exhibits many 
of the features found from molecular dynamic simulations |im ITT] . Some new 
features not previously described have also been found and it should be pos- 
sible to check if these features also exist in conventional molecular dynamic 
simulations. 



2 The Hydrodynamic Equations 

In this paper, we used the hydrodynamic equations postulated by Haff|l]. 
They are essentially the Navier-Stokes equations with the transport coeffi- 
cients being dependent upon the density and the granular temperature [2] 
and with an extra dissipation term due to the inelastic nature of the particles. 
Written out fully, the equations are[H] 

dp 

m 

dt 
dT 

'dt 



where sums over repeated indices is implied, p stands for the density, Ui for 
the components of the velocity, t for the time, T for the granular temperature |13j. 
Pf for the pressure due to the free energy, rj for the viscosity tensor, k, for 
the thermal conductivity and 7 for the thermal dissipation coefficient. The 
derivation of these equations is discussed in detail in Haff |1] ^1] . We use the 
usual form 

Vijki = vi^ikSji + Sii6jk) + x^ij^ki (4) 
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for the viscosity tensor where t] and x stand for the shear and bulk viscosities 
respectively. 

Similar equations can be derived JUIIT^ by performing a Chapman- Enskog 
expansion of the Boltzmann equations for an inelastic hard sphere model 
about the homogeneous state. These equations are rather similar except for 
some scaling differences for 7 and rj and an extra term involving the density 
gradient. The pressure term in this model is just of the ideal gas form. The 
viscosity also has a very simple form with the bulk viscosity being equal to 
the shear viscosity. 

Following Hill and Mazenko|H], we simulated the Haff equations with a 
pressure term which reflects the interactions between the grains. 

3 The simulation of the hydrodynamic equa- 
tions 

The Haff equations were simulated with the same basic algorithm as in Hill 
and Mazenko^. The algorithm was improved with the use of a superior 
adaptive step technique which speeded up the simulations significantly with- 
out loss of accuracy. As in Hill and Mazenko |8j , the pressure term is chosen 
in the model to account for the excluded volume effects in the granular fluid. 
The pressure, is, as usual, related to the free energy density / by p/ = p^ — /. 
The interaction part of the free energy density used for this purpose was cho- 
sen to be of the form 

f{p) = K{p-p,re{p~p,) (5) 

where is the step function and pc is the packing density. In Hill and 
MazenkofS], a = 2. This free energy corresponds to particles which behave 
roughly like stiff springs with regard to compression. In our work, we chose 
a = 10 to preserve the continuity of the higher derivatives of f{p) which 
enables the higher order accuracy of the fourth-order Runge-Kutta method 
to be exploited. This roughly corresponds to particles with a very thin soft 
exterior and a hard interior. In the simulations, it is however not difficult to 
verify that the precise form of this free energy density is not important as 
long as the interaction is stiff enough. The large value of a permits a much 
smaller range of densities above the packing density. The value K = A x 10^ 
was chosen in Hill and MazenkojH] while we chose K = 4x 10^^ to adjust for 
the larger value of a. 
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Figure 1: The x component of the velocity for a system that has been at- 
tracted to a final shear state (r^ = 5, 7 = 1, t = 10^). 



We also added a fictional term to the free energy density 



(6) 



where p^ = 10^'^, A = 10^°. This contribution to the pressure prevents areas 
of extremely low density from being produced since these cause numerical 
instabilities. 



4 Observations 

The units used in this paper are the same as in Ref. [5]. The initial con- 
ditions chosen for most of the simulations were a uniform density half the 
packing density {pc = 0.2), random uniformly distributed initial velocities 
in the range [—1, 1] for each component (ensuring that the initial center of 
mass velocity was zero) and a uniform granular temperature of 1. These 
conditions were chosen for their relative simplicity. The precise random con- 
figuration of velocities was seen in some cases to matter significantly in the 
subsequent dynamical evolution of the system including the form of the final 
state achieved. The dynamical equations were solved numerically on a grid 
in space where the lattice spacing was chosen to be 0.9 for most simulations, 
a value which was verified to provide good accuracy for the values of viscos- 
ity that we used. The densities were stored at the intersection points of the 




Figure 2: The x and z velocity distributions for a vortex state seen as a 
greyscale map (intensity of white proportional to absolute velocity, 77 = 9, 7 = 
11, t = 10^) . For the x velocity distribution shown in the left image, the 
large white region on the right has a negative velocity while for the z velocity 
distribution shown in the right, the large white region on the left has a 
positive velocity. 



grid while the velocities were stored at the mid-points of the lines between 
the intersections (x-components on the lines perpendicular to the x-axis and 
the z-components on the lines perpendicular to the z-axis). The upwind 
technique was then used to determine which velocity corresponded to which 
density. The method is described in detail in Ref. j^. The fourth order 
Runge-Kutta method was used for the time evolution and an adaptive time 
step based on the estimated error from a single time step and two steps of 
half the size used to control the error and speed up the simulation. 

There are a number of parameters in the model. It is important to realize 
that the viscosity sets the major length scale in the problem and larger 
viscosities enable us to use larger lattice spacings without introducing large 
numerical errors. For some low viscosity simulations, we have used smaller 
lattice spacings. 

In contrast to the work of McNamara and Young[T7j. our simulations have 
been run long enough to access the final attracting states. The simulations 
have also been performed for a wide range of values in parameter space. 
These final attracting states have been discussed in some detail in Soto and 
Mareschal|18j. There is the well known homogenous cooling state (HCS) 
with uniformly decreasing average velocity and temperature, and no structure 
formation. There is the shear band configuration^] with the velocity field 
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Figure 3: The x and z velocity distributions for a system in the vortex state 
(r^ = 9,7 = ll,t = 109). 



taking on the form v = Vo{t) cos(27ra;/L)z or v = Vo{t) cos{2TczL)St where z 
and X are the unit vectors in the z and x directions respectively and L is the 
size of the domain. An example of the distribution of the x component of 
the velocity distribution in a shear band state in the x direction is shown in 
figure HJ 

The final type of attractor is a cyclical one with features similar to a 
vortex structure. If we define 
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m 



(7) 



where E(t) = {pu^ + T) stands for the energy at time t and consider the 
rescaled velocities v = u/'y and differentially rescaled time 



ds = 'ydt 



(8) 



we find that f is a periodic function of s. A greyscale map of the x and 
z components of the velocity distribution at a point in the cycle is shown 
in figure El Seen this way, the cyclical structure looks like two tadpole-like 
structures moving much like two shear bands. This seems to be related to 
the state with two counter-rotating vortices found in the nonlinear analysis of 
Soto and Mareschal_18j where the velocity is dominated by the kx = ^,kz = 
^ mode in Fourier space. Another view of the velocity distribution is s hown 
in figure El The density distribution at various times in the cycle is shown 
in figure m 
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Figure 4: The density distribution for the vortex state (density is propor- 
tional to the intensity of white in the figure). The time goes from left to 
right and top to bottom, with the values being t = 1.734 x 10*^, 1.842 x 
10^ 2.197 X 10^ 2.895 x 10^ 3.779 x 10^ 4.545 x 10^ 5.437 x 10^ 7.252 x 10^ 
and 9.391 x 10^ 

We also find, as in past studies [T71 ITUl ITT] , an intermediate clustering 
state which is metastable. However, which particular final attracting state is 
chosen seems to depend on many criteria including the starting distribution 
of velocities. The homogeneous cooling state is chosen by the system when 
the dissipation is low|18j. However, if the system size is large enough, any 
dissipation will make the homogeneous cooling state unstable jTHj. Soto and 
Mareschal[18^ carried out a non- linear stability analysis for the shearing and 
vortex states and found a criterion for the stability of the states provided 
the deviations from the homogeneous cooling state were small. The criterion 
is relatively simple but is accurate only when the state does not deviate too 
much from the homogeneous cooling state. In our simulations, shear bands 
were mostly observed but some vortex states were also observed even with 
the same parameters but different system sizes or initial conditions. For 
simulations with k = 10, x = 1, the vortex attracting state was achieved for 
(77,7) = (7, 7), (7, 15), (9, 11) and (9, 19) in one of the runs (since the final 
state does depend on the initial conditions, it occured for different values of 
the parameters in other runs). It is also interesting to note that the velocity 
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Figure 5: A plot of logE versus logi for a system that never leaves the HCS 
[t] = 25,7 = 0.1). A line with slope -2 is drawn for reference purposes. 

field is always the first to go to the attracting state with the other fields 
following much later. 

We can be somewhat more quantitative if we look at the evolution of 
the energy (defined as E = {p{\u'^ + ^))) ^^^ the granular thermal energy 
(defined as E^ = (pT)) with time|171 ITTH ITT] . It is well known that in the 
HCS, the energy evolves as 

so that at large times, the slope of a log-log plot of energy with respect to 
time has a slope of -2. We show this plot for a system which never leaves the 
HCS in figure El The plot for a system which has the clustering instability 
is shown in figure El We see that during the clustering instability, the rate 
of decrease of energy is reduced. Some studies jTT] indicate that the rate of 
decrease of the energy goes as -E ~ t^^ during this metastable phase but we 
find that the relation depends on the value of the parameters. The plot for a 
system which has the clustering instability and which settles into the vortex 
state is shown in figure [71 The reduction in the rate of decrease of energy 
in the clustering phase is similar to that of a system that finally goes into a 
shear band state but the cyclical nature of the final state is clearly seen in 
this plot. The average slope of the plot in the final state is still -2. 

The evolution of the granular thermal energy Ej- = (pT) has a similar 
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Figure 6: A plot of logE versus logt for a system that has a clustering 
instability and which finally settles into a shear band state (r^ = 5, 7 = 1). A 
line with slope -2 is drawn for reference purposes. 

behaviour. The corresponding log-log plots are shown in figures 151 lUl and ITUl 
It is very interesting to note that the granular thermal energy Ej- does go as 
t^^ for the initial evolution of the clustering state and this does seem to be 
independent of the parameters. This result might be related to the one found 
by Xiaobo et alfTT^ and Miller and Luding 19, where the kinetic energy decays 
as t~^ in molecular dynamics simulations. The slope of the plots, discussed 
in some detail in Hill and MazenkojEj, Ve = ^^ j^ , is another interesting 
quantity. Examples of the behaviour of Ve are shown for several values of 
the parameters in figures ^2 ^1 El and El respectively. It is interesting 
to note that in systems where there is a clustering instability, Ve has large 
fluctuations in the later part of the evolution of this phase. 

Since the energy behaves as Q in the HCS and the system is initially in 
the HCS, we have 

Since to is relatively small (of order 1), Ve should rapidly decrease from to 
-2 if the system stays in the HCS. Any increase in Ve signals a departure from 
the HCS and we can hence identify the time at which Ve reaches a minimum 
as the time the system begins to depart from the HCS. We define this time 
to be tE- It was found earlier jH] that tE^^ followed a power law with respect 
to 7, that is tE ~ 7'^. This result cannot hold for all 7 since, for a given 



Figure 7: A plot of logE versus logt for a system that has a clustering 
instability and which finally settles into a vortex state (77 = 9,7 = 11). 

system size, there is a finite cuttoff 7c below which the system never leaves 
the HCS. However, for 7 well above this value (which is generally much less 
than 1 except for high viscosities), we find a power law relation for several 
values of r] (x = 1, K = 10 being held fixed) with the exponent /i being in the 
range [—1.5, —1.4]. Examples of the power law behavior for several values of 
the viscosities are shown in figure IT^ . 

Another quantity of interest which has been analysed earlier in the literature jTT] 
is the ratio of the thermal energy to the kinetic energy Et/ Ek. For a vortex 
state, the behaviour is interesting in that this ratio oscillates about a fixed 
value as seen in figure E| This is another example of the cyclical behaviour 
of this attractor. For a shear state, the ratio approaches a constant value 
asymptotically as seen in figure IT71 

We now look at the clustering instability a bit more closely so as to 
identify clearly when it occurs. Since phase separation can be due to the 
clustering instability or a result of the velocity gradients in an attracting 
state, we have to distinguish between these two. To do so, we characterize 
phase separation as having occured if the rms fluctuation in the density is 
above 0.3 and the system has reached the packing density at at least one 
point. We call the time at which this first occures as tphase- We then identify 
the time t/maz at which the system has reached its final state. The main 
feature of the shear band state is the picking of one direction in which the 
bands occur. The result of this is a dominance of kinetic energy in that 
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Figure 8: A plot of log £"7 versus logt for a system that never leaves the HCS 
(?7 = 25,7 = 0.1). A line with slope -2 is drawn for reference purposes. 

direction. Hence, for a shear band state, we identify t final as the first time at 
which the kinetic energy in one direction dominates the kinetic energy in the 
other direction by a factor of 1000, that is the first time when j^^ > 1000 

or when -(—^ ^ 1000. A similar criterion is used to identify the shearing 
state in McNamara and Young [T7j. For a vortex state, we indentify t final as 
the time at which the oscillations in Er/Ek start. If tphase < tfinai we claim 
that the phase separation is a result of the clustering instability. Using this 
criterion, the phase diagram for the clustering instability is shown in figure 



5 Summary and Conclusions 

To summarize, we have used a hydrodynamic model of granular dynamics to 
simulate the cooling of a granular gas with no external field and with periodic 
boundary conditions over a large region of parameter space. We have found 
that the results of simulations using this model are in broad agreement with 
molecular dynamic simulations of the same phenomenon. We find that the 
system is attracted to one of three final states - the homogeneous cooling 
state, the shear band state or a vortex state. The HCS is the final state only 
if the system is small and the dissipation is very low while the shear band 
or vortex states result when the dissipation is high with the exact criterion 
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Figure 9: A plot of logi^r versus logt for a system that has a clustering 
instability and which finally settles into a shear band state (77 = 5,7 = 1). 
Lines with slope -1 and -2 are drawn for reference purposes. 

for the selection of the state unknown though it has been found to depend 
on the initial state of the system. We have also found a phase diagram for 
the clustering instability in this model and have confirmed the result in Hill 
and Mazenko|H] that the time taken to depart from the homogeneous cooling 
state follows a power law with respect to the dissipation coefficient (for high 
enough dissipation) over a much larger region of parameter space. It should 
be possible to check this result with molecular dynamic simulations. 
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Figure 12: Ve versus time for r] = 25, x = 1, 7 = 19, /t = 10. This state goes 
finally to the shear state but only after a clustering instability phase. 




Figure 13: Ve versus time for r^ = 7, x = 1, 7 = 19, k = 10. This state goes 
finally to the shear state but only after a strong clustering instability phase. 
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Figure 14: The derivative of the logarithm of the thermal energy with respect 
to the logarithm of the time for a system with rj = 25, x = 1, /« = 10, 7 = 0.1 
which does not exit the HCS. We see that the derivative settles to -2 and 
does not change. The small deviations from -2 seem to be due to numerical 
errors. 
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Figure 15: The power law Ie = 7^* for rj = 5,7,11,25. Note that for r/ = 
11, 25, we have to exclude the first point from the fit since it is too close to 
7c. The values for /x are -1.39, -1.50, -1.41 and -1.40 respectively. 
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Figure 16: The ratio of the thermal to kinetic energy Et/ E^ for r^ = 9, 7 = 
11, X = I,/? = 10 which goes to a final vortex state. Note that once the 
system reaches the vortex state, the ratio oscillates about a constant value. 




Figure 17: The ratio of the thermal to kinetic energy Et/ Ek for 77 = 5,7 = 
3, X = 1, /t = 10 which goes to a final shear state. Note that once the system 
reaches the shear state, the ratio converges to a constant. 
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Figure 18: The phase diagram for the clustering instabihty. The values of 
the other parameters are x = 1, 't = 10. Clustering occurs for parameters to 
the right of the dashed line. 
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Rescaled x component of velocity for eta=7, chi=1 , kappa=1 
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Rescaled z component of velocity for eta=7, chi=1 , kappa=1 
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